# This analysis was conducted using R Studio Version 1.1.419 for Mac and R version 3.5.2 for Mac
sink("./logofcode.txt")
rm(list = ls(all = TRUE))
install.packages(pacman)
pacman::p_load(foreign, ggplot2, plm, reshape2, countrycode, sandwich, lmtest, MASS, 
               rworldmap, RColorBrewer, states, mice, VIM, stargazer, margins, clusterSEs, lme4, optimx,
               coefplot, lattice, survey, dplyr, wordcloud, sampleSelection, quanteda,
               Matrix, ldatuning, topicmodels, readtext, stm, lda, bursts, tidytext, 
               tm, readstata13, fastDummies, panelView, igraph, NetMix, dnr, statnet, network, sna, 
               DataCombine, sensemakr, estimatr, grid, gridExtra)
setwd("ENTER DIRECTORY") # enter your directory

####################### Create Datasets (skipped for replication log) ####
##->## Main dataset no lag
IMF <- read.csv("IMF_cond.csv", fileEncoding="UTF-8-BOM") # Load the conditions tab from Kentikelenis et al. (2016) downloaded in October 2018
IMF$count <- 1 # create a counter variable equal to 1 for each condition
IMF$Waiver[is.na(IMF$Waiver)] <- 0 # Recode no waiver as 0
IMF$Waiver[IMF$Waiver > 0] <- 1 # Recode all waivers as 1
for (i in 1:32261){
  IMF$DEB[i] <- ifelse(IMF$Policy.Area..short.[i] == "DEB", 1, 0)
  IMF$FIN[i] <- ifelse(IMF$Policy.Area..short.[i] == "FIN", 1, 0)
  IMF$FP[i] <- ifelse(IMF$Policy.Area..short.[i] == "FP", 1, 0)
  IMF$EXT[i] <- ifelse(IMF$Policy.Area..short.[i] == "EXT", 1, 0)
  IMF$RTP[i] <- ifelse(IMF$Policy.Area..short.[i] == "RTP", 1, 0)
  IMF$SOE[i] <- ifelse(IMF$Policy.Area..short.[i] == "SOE", 1, 0)
  IMF$LAB[i] <- ifelse(IMF$Policy.Area..short.[i] == "LAB", 1, 0)
  IMF$PRI[i] <- ifelse(IMF$Policy.Area..short.[i] == "PRI", 1, 0)
  IMF$SP[i] <- ifelse(IMF$Policy.Area..short.[i] == "SP", 1, 0)
  IMF$POV[i] <- ifelse(IMF$Policy.Area..short.[i] == "POV", 1, 0)
  IMF$INS[i] <- ifelse(IMF$Policy.Area..short.[i] == "INS", 1, 0)
  IMF$ENV[i] <- ifelse(IMF$Policy.Area..short.[i] == "ENV", 1, 0)
  IMF$OTH[i] <- ifelse(IMF$Policy.Area..short.[i] == "OTH", 1, 0)
} # code each condition by policy category
IMFagg <- aggregate(list(IMF$count, IMF$DEB, IMF$FIN, IMF$FP, IMF$EXT, IMF$RTP, IMF$SOE, IMF$LAB, 
                         IMF$PRI, IMF$SP, IMF$POV, IMF$INS, IMF$ENV, IMF$OTH, IMF$Waiver), by = list(IMF$Country.Name, 
                                                                                         IMF$Country.Code, 
                                                                                         IMF$Arrangement.Date, 
                                                                                         IMF$Arrangement.ID, 
                                                                                         IMF$Year, 
                                                                                         IMF$Arrangement.Duration,
                                                                                         IMF$Arrangement.Type,
                                                                                         IMF$Arrangement.Amount..SDR), 
                    FUN = sum) # aggregate conditions by category and summing by project

colnames(IMFagg) <- c("ctry", "ccode", "date", "ID", "year", "dur", "type", "loansize", "count", "DEB", "FIN", "FP", "EXT", "RTP", 
                      "SOE", "LAB", "PRI", "SP", "POV", "INS", "ENV", "OTH", "waiver") # change column names
IMFagg$dur <- gsub("(.*),.*", "\\1", IMFagg$dur) #Eliminate anything to right of comma in duration column
IMFagg$type <- gsub("(.*),.*", "\\1", IMFagg$type) #Eliminate anything to right of comma in type column
IMFagg$type <- gsub("(.*),.*", "\\1", IMFagg$type) #Repeat for those with two commas
IMFagg$PRGF <- 0 # create blank PRGF column
for (i in 1:2121) {
  IMFagg$PRGF[i] <- ifelse(IMFagg$type[i]  == "PRGF"  | IMFagg$type[i] == "ESAF" | 
                           IMFagg$type[i]  == "SAF" | IMFagg$type[i]  == "ECF" | 
                           IMFagg$type[i]  == "RCF" | IMFagg$type[i]  == "SCF", 1, 0)
} # code PRGF observations
write.csv(IMFagg, "IMF_fin.csv") # write the full conditionality data to file

WDI <- read.csv("WDI10.24.19.csv", fileEncoding="UTF-8-BOM") # Load data for economic covariates downloaded from World Development Indicators in October 2019
WDI <- WDI[-c(4,13)] # Delete unnecessary columns
colnames(WDI) <- c("ctry", "ccode", "year", "gdp", "gdppc", "USaid", "DACaid", "trade", "dservexp", "dshortexp", "reserves", "curraccgdp",
                   "fdigdp", "inflation") # rename columns
WDI[WDI == ".."] <- NA # change two dot missing value notation to NA
WDI$USaid <- as.numeric(as.character(WDI$USaid)) # change to numeric
WDI$USaid <- WDI$USaid / 1000000 # convert to millions
WDI$USaid[WDI$USaid < 0] <- 0 # change negatives to 0
WDI$DACaid <- as.numeric(as.character(WDI$DACaid)) # change to numeric
WDI$DACaid <- WDI$DACaid / 1000000 # convert to millions
WDI$DACaid[WDI$DACaid < 0] <- 0 # change negatives to 0
WDI <- WDI[-1] # delete country name (would be duplicate)
prelim <- merge(IMFagg, WDI, by = c("ccode", "year"), all.x = T) # merge data

quota <- read.csv("quota.csv", fileEncoding="UTF-8-BOM") # load data on quotas downloaded from IMF website in October 2018
colnames(quota) <- c("ctry", 1974:2014) # change column names
quota <- melt(quota, id.vars=c("ctry")) # switch to column format instead of rows
colnames(quota) <- c("ctry", "year", "quota") # change column names
quota$ccode <- countrycode(quota$ctry, origin = "country.name", destination = "iso3c") # build country code
quota$quota <- quota$quota / 1000000 # convert to millions
quota <- quota[-1] # delete country name (would be duplicate)
prelim <- merge(prelim, quota, by = c("ccode", "year"), all.x = T) # merge data

load("Dyadicdata.RData") # load dyadic UN voting data downloaded from Bailey et al. (2017) in November 2018 V21 https://doi.org/10.7910/DVN/LEJUQZ
UNGA <- x # assign loaded data to name
UNGA_US <- UNGA[UNGA$ccode1 == 2,] # subset to relations with US
UNGA_US$ccode <- countrycode(UNGA_US$ccode2, origin = "cown", destination = "iso3c") # build ccode
UNGA_US <- UNGA_US[-c(1:2,4:5,7:21)] # delete excess columns
UNGA_US <- aggregate(list(UNGA_US$absidealdiff, UNGA_US$absidealimportantdiff), 
                     by = list(UNGA_US$ccode, UNGA_US$year), FUN = mean) # eliminates issue with duplicates
colnames(UNGA_US) <- c("ccode", "year", "absidealdiff", "absidealimportantdiff")
prelim <- merge(prelim, UNGA_US, by = c("ccode","year"), all.x = T) # merge data

polity <- read.csv("p4v2017.csv", fileEncoding="UTF-8-BOM") # load polity2 v 2017 data downloaded from systemic peace in November 2018
# note that I deleted extra columns in Excel, leaving only the columns named "ccode", "country", "year", and "polity2"
# which appear as column numbers 2, 4, 5, and 11 in the original dataset. "Country" was renamed "ctry" to match my other codings
polity$ccode <- countrycode(polity$ctry, origin = "country.name", destination = "iso3c") # build ccode
polity <- polity[-c(2)] # delete unnecessary column
polity <- aggregate(list(polity$polity2), 
                     by = list(polity$ccode, polity$year), FUN = mean) # eliminates issue with duplicates
colnames(polity) <- c("ccode", "year", "polity2") # change column names
prelim <- merge(prelim, unique(polity), by = c("ccode","year"), all.x = T) # merge data

dpi <- read.csv("DPI.csv", fileEncoding="UTF-8-BOM") # load DPI data on number of checks downloaded in November 2018 https://publications.iadb.org/en/database-political-institutions-2017-dpi2017
dpi$ccode <- countrycode(dpi$countryname, origin = "country.name", destination = "iso3c") # add country code
dpi <- dpi[-c(1)] # drop country name column
dpi <- aggregate(list(dpi$checks), by = list(dpi$ccode, dpi$year), FUN = max) # eliminates issue with duplicates
colnames(dpi) <- c("ccode", "year", "checks") # change column names
prelim <- merge(prelim, unique(dpi), by = c("ccode", "year"), all.x = T) # merge data

last <- read.csv("IMFlast.csv", fileEncoding="UTF-8-BOM") # load hand-coded data on time since last IMF program; I utilized the 
# Kentikelenis et al. (2016) conditionality data imported above to manually calculate the number of 
# years since a country's last IMF program data
last <- aggregate(list(last$last), by = list(last$ccode, last$year), FUN = max) # eliminates issue with duplicates
colnames(last) <- c("ccode", "year", "last") # change column names
prelim <- merge(prelim, unique(last), by = c("ccode", "year"), all.x = T) # merge data

nelson <- read.dta("nelson.dta") # load replication data from Nelson (2014) on liberal ideology among IMF recipients https://doi-org.ezproxy.cul.columbia.edu/10.1017/S0020818313000477
nelson <- nelson[-c(2,4:29, 31:49)] # eliminate extra columns
colnames(nelson) <- c("ccode", "year", "prop_neoliberal") # rename columns
nelson$ccode <- countrycode(nelson$ccode, origin = "cown", destination = "iso3c") # harmonize ccodes
prelim <- merge(prelim, unique(nelson), by = c("ccode", "year"), all.x = T) # merge data

rfa <- read.csv("RFA_ctry.csv", fileEncoding="UTF-8-BOM") # load hand coded data on membership in and contributions from cooperative/competitive RFAs; coding is discussed in the paper's Appendix 3
rfa$ccode <- countrycode(rfa$ctry, origin = "country.name", destination = "iso3c") # generate ccodes
rfa <- rfa[-1] # eliminate duplicate ctry column
prelim <- merge(prelim, unique(rfa), by = c("ccode", "year"), all.x = T) # merge
prelim$ooamt[is.na(prelim$ooamt)] <- 0 # set NAs which represent no membership to 0 
prelim$oomem[is.na(prelim$oomem)] <- 0 # set NAs which represent no membership to 0 
prelim$cofin[is.na(prelim$cofin)] <- 0 # set NAs which represent no membership to 0 
prelim$infoshare[is.na(prelim$infoshare)] <- 0 # set NAs which represent no membership to 0 
prelim$coopmem[is.na(prelim$coopmem)] <- 0 # set NAs which represent no contribution to 0 
prelim$coopamt[is.na(prelim$coopamt)] <- 0 # set NAs which represent no contribution to 0 
prelim$AMFm[is.na(prelim$AMFm)] <- 0 # set NAs which represent no membership to 0 
prelim$CRAm[is.na(prelim$CRAm)] <- 0 # set NAs which represent no membership to 0 
prelim$CMIMm[is.na(prelim$CMIMm)] <- 0 # set NAs which represent no membership to 0 
prelim$EFSDm[is.na(prelim$EFSDm)] <- 0 # set NAs which represent no membership to 0 
prelim$EFSMm[is.na(prelim$EFSMm)] <- 0 # set NAs which represent no membership to 0 
prelim$EUBOPm[is.na(prelim$EUBOPm)] <- 0 # set NAs which represent no membership to 0 
prelim$ESMm[is.na(prelim$ESMm)] <- 0 # set NAs which represent no membership to 0 
prelim$FLARm[is.na(prelim$FLARm)] <- 0 # set NAs which represent no membership to 0 

nelda <- read.csv("NELDA.csv", fileEncoding="UTF-8-BOM") # load election data from NELDA in November 2018
for (i in 1:3104){
  nelda$leg[i] <- ifelse(nelda$types[i] == "Legislative/Parliamentary", 1, 0)
  nelda$pres[i] <- ifelse(nelda$types[i] == "Executive", 1, 0)
} # code elections by type
nelda <- nelda[-c(2:4,6:125)] # delete unnecessary rows
colnames(nelda)[1] <- "ccode" # change column name
nelda <- aggregate(list(nelda$leg, nelda$pres), by = list(nelda$ccode, nelda$year), FUN = sum) # aggregate by election type
colnames(nelda) <- c("ccode", "year", "leg", "pres") # change column names
nelda$leg[nelda$leg > 1] <- 1 # change to binary
nelda$pres[nelda$pres > 1] <- 1 # change to binary
nelda$elec <- nelda$leg + nelda$pres # create overall election measure by summing
nelda$elec[nelda$elec > 1] <- 1 # change to binary
prelim <- merge(prelim, unique(nelda), by = c("ccode", "year"), all.x = TRUE) # merge data
prelim$elec[is.na(prelim$elec)] <- 0 # change NA to 0 because no election is currently set as blank
prelim$pres[is.na(prelim$pres)] <- 0 # change NA to 0 because no election is currently set as blank
prelim$leg[is.na(prelim$leg)] <- 0 # change NA to 0 because no election is currently set as blank

ucdp <- read.csv("ucdp.csv", fileEncoding="UTF-8-BOM") # load conflict/war data downloaded from UCDP in November 2018
ucdp <- ucdp[-c(4:18)] # Delete unnecessary columns
states <- gwstates # load Gleditch and Ward list of independent states (included in R)
states$ccode <- countrycode(states$iso3c, origin = "cowc", destination = "iso3c") # generate country code
states <- states[-c(2:6)] # Delete unnecessary columns
colnames(states) <- c("gwno", "ccode") # Change column names
ucdp2 <- merge(ucdp, states, by = c("gwno")) # merge data
ucdp2 <- ucdp2[-c(1,4)] # delete extra columns
colnames(ucdp2) <- c("year", "war", "ccode") # change column names
prelim <- merge(prelim, unique(ucdp2), by = c("ccode", "year"), all.x = TRUE) # merge data

unsc <- read.csv("unsc.csv", fileEncoding="UTF-8-BOM") # load data on UNSC membership from Dreher, Sturm, and Vreeland (2009) 
# downloaded in November 2018; I hand-coded the more recent years from UN online materials to bring the data up to date
unsc$ccode <- countrycode(unsc$ctry, origin = "country.name", destination = "iso3c") # harmonize country codes
unsc$unsc[unsc$unsc == "."] <- 0 # change dot notation for missing values to NA
unsc <- na.omit(unsc) # eliminates countries like East Germany
unsc <- unsc[-c(1,4)] # deletes unnecessary columns
prelim <- merge(prelim, unique(unsc), by = c("ccode", "year"), all.x="T") # merge data
prelim$unsc[is.na(prelim$unsc)] <- 0 # NAs are 0s
write.csv(prelim, "prelim.csv") # write full data to csv

##->## Main data no lag with imputation
prelim$ccode <- as.character(prelim$ccode) # change various data formats for imputation purposes
prelim$gdp <- as.numeric(as.character(prelim$gdp))
prelim$gdppc <- as.numeric(as.character(prelim$gdppc))
prelim$trade <- as.numeric(as.character(prelim$trade))
prelim$dservexp <- as.numeric(as.character(prelim$dservexp))
prelim$dshortexp <- as.numeric(as.character(prelim$dshortexp))
prelim$reserves <- as.numeric(as.character(prelim$reserves))
prelim$curraccgdp <- as.numeric(as.character(prelim$curraccgdp))
prelim$fdigdp <- as.numeric(as.character(prelim$fdigdp))
prelim$inflation <- as.numeric(as.character(prelim$inflation))
prelim$unsc <- as.numeric(as.character(prelim$unsc))
prelim$loansize <- as.numeric(as.character(prelim$loansize))
prelim$dur <- as.numeric(prelim$dur)
prelimimp <- mice(prelim, method = "cart", seed = 500, m = 1,
                  pred=quickpred(prelim, exclude=c("X", "ID", "ctry", "date", "type"))) # impute missing data
prelimimp <- complete(prelimimp) # extract imputed values to data frame
write.csv(prelimimp, "prelimimp.csv") # write full data to csv

##->## Selection data with lags and imputation
mems <- read.csv("members.csv", fileEncoding="UTF-8-BOM") # load hand coded data on IMF and RFA memberships over time
mems$ccode <- countrycode(mems$ctry, origin = "country.name", destination = "iso3c") # generate country codes
part <- data.frame(cbind(as.character(prelim$ccode), as.character(prelim$year))) # create data frame with active project country-years
part$part <- 1 # create participation variable
colnames(part) <- c("ccode", "year", "part") # change column names
part$year <- as.numeric(as.character(part$year)) # change to numeric
part$ccode <- as.character(part$ccode) # change to character
mems <- merge(mems, part, by = c("ccode", "year"), all.x=T) # merge data
mems$part[is.na(mems$part)] <- 0 # change non-participating rows to 0
WDI$year <- WDI$year + 1 # lag by one year
mems <- merge(mems, WDI, by = c("ccode", "year"), all.x = T) # from here down, I am just repeating the above steps for dataset creation with the selection data
quota$year <- as.numeric(as.character(quota$year)) + 1 # lag by one year
mems <- merge(unique(mems), unique(quota), by = c("ccode", "year"), all.x = T)
UNGA_US$year <- UNGA_US$year + 1
mems <- merge(mems, UNGA_US, by = c("ccode", "year"), all.x = T)
polity$year <- polity$year + 1
mems <- merge(mems, unique(polity), by = c("ccode", "year"), all.x = T)
last$year <- last$year + 1
mems <- merge(mems, last, by = c("ccode", "year"), all.x = T)
rfa$year <- rfa$year + 1
mems <- merge(mems, unique(rfa), by = c("ccode", "year"), all.x = T)
mems$ooamt[is.na(mems$ooamt)] <- 0
mems$coopamt[is.na(mems$coopamt)] <- 0
mems$oomem[is.na(mems$oomem)] <- 0 
mems$coopmem[is.na(mems$coopmem)] <- 0 
mems$cofin[is.na(mems$cofin)] <- 0 
mems$infoshare[is.na(mems$infoshare)] <- 0 
mems$AMFm[is.na(mems$AMFm)] <- 0 
mems$CRAm[is.na(mems$CRAm)] <- 0 
mems$CMIMm[is.na(mems$CMIMm)] <- 0 
mems$EFSDm[is.na(mems$EFSDm)] <- 0 
mems$EFSMm[is.na(mems$EFSMm)] <- 0
mems$EUBOPm[is.na(mems$EUBOPm)] <- 0
mems$ESMm[is.na(mems$ESMm)] <- 0
mems$FLARm[is.na(mems$FLARm)] <- 0 
nelda$year <- nelda$year + 1
mems <- merge(mems, nelda, by = c("ccode", "year"), all.x=T)
mems$elec[is.na(mems$elec)] <- 0 
mems$pres[is.na(mems$pres)] <- 0 
mems$leg[is.na(mems$leg)] <- 0
ucdp2$year <- ucdp2$year + 1
mems <- merge(mems, unique(ucdp2), by = c("ccode", "year"), all.x=T)
unsc$year <- unsc$year + 1
mems <- merge(mems, unique(unsc), by = c("ccode", "year"), all.x=T)
mems$unsc[is.na(mems$unsc)] <- 0

mems$ccode <- as.character(mems$ccode)
mems$gdp <- as.numeric(as.character(mems$gdp))
mems$gdppc <- as.numeric(as.character(mems$gdppc))
mems$trade <- as.numeric(as.character(mems$trade))
mems$dservexp <- as.numeric(as.character(mems$dservexp))
mems$dshortexp <- as.numeric(as.character(mems$dshortexp))
mems$reserves <- as.numeric(as.character(mems$reserves))
mems$curraccgdp <- as.numeric(as.character(mems$curraccgdp))
mems$fdigdp <- as.numeric(as.character(mems$fdigdp))
mems$inflation <- as.numeric(as.character(mems$inflation))
mems$unsc <- as.numeric(as.character(mems$unsc))
memsimp <- mice(mems, method = "cart", seed = 500, m = 1,
                  pred=quickpred(mems, exclude=c("X", "ctry")))
memsimp <- complete(memsimp)
write.csv(memsimp, "memsimp.csv")

##<-## Main dataset with lag (some coviariates still lagged from above, so coding is adjusted below)
prelimlag <- IMFagg # steps are again the same as above
prelimlag <- merge(prelimlag, WDI, by = c("ccode", "year"), all.x = T)
prelimlag <- merge(prelimlag, quota, by = c("ccode", "year"), all.x = T) 
prelimlag <- merge(prelimlag, UNGA_US, by = c("ccode","year"), all.x = T) 
prelimlag <- merge(prelimlag, unique(polity), by = c("ccode","year"), all.x = T) 
dpi$year <- dpi$year + 1
prelimlag <- merge(prelimlag, unique(dpi), by = c("ccode", "year"), all.x = T)
prelimlag <- merge(prelimlag, unique(last), by = c("ccode", "year"), all.x = T)
nelson$year <- nelson$year + 1
prelimlag <- merge(prelimlag, unique(nelson), by = c("ccode", "year"), all.x = T)
prelimlag <- merge(prelimlag, unique(rfa), by = c("ccode", "year"), all.x = T)
prelimlag$ooamt[is.na(prelimlag$ooamt)] <- 0 
prelimlag$oomem[is.na(prelimlag$oomem)] <- 0
prelimlag$coopmem[is.na(prelimlag$coopmem)] <- 0 
prelimlag$coopamt[is.na(prelimlag$coopamt)] <- 0 
prelimlag$cofin[is.na(prelimlag$cofin)] <- 0 
prelimlag$infoshare[is.na(prelimlag$infoshare)] <- 0 
prelimlag$AMFm[is.na(prelimlag$AMFm)] <- 0 
prelimlag$CRAm[is.na(prelimlag$CRAm)] <- 0
prelimlag$CMIMm[is.na(prelimlag$CMIMm)] <- 0 
prelimlag$EFSDm[is.na(prelimlag$EFSDm)] <- 0  
prelimlag$EFSMm[is.na(prelimlag$EFSMm)] <- 0 
prelimlag$EUBOPm[is.na(prelimlag$EUBOPm)] <- 0 
prelimlag$ESMm[is.na(prelimlag$ESMm)] <- 0 
prelimlag$FLARm[is.na(prelimlag$FLARm)] <- 0 
prelimlag <- merge(prelimlag, unique(nelda), by = c("ccode", "year"), all.x = TRUE)
prelimlag$elec[is.na(prelimlag$elec)] <- 0
prelimlag$pres[is.na(prelimlag$pres)] <- 0 
prelimlag$leg[is.na(prelimlag$leg)] <- 0
prelimlag <- merge(prelimlag, unique(ucdp2), by = c("ccode", "year"), all.x = TRUE)
prelimlag <- merge(prelimlag, unique(unsc), by = c("ccode", "year"), all.x="T")
prelimlag$unsc[is.na(prelimlag$unsc)] <- 0
write.csv(prelimlag, "prelimlag.csv")

##->## Main data with lag and imputation
prelimlag$ccode <- as.character(prelimlag$ccode)
prelimlag$gdp <- as.numeric(as.character(prelimlag$gdp))
prelimlag$gdppc <- as.numeric(as.character(prelimlag$gdppc))
prelimlag$trade <- as.numeric(as.character(prelimlag$trade))
prelimlag$dservexp <- as.numeric(as.character(prelimlag$dservexp))
prelimlag$dshortexp <- as.numeric(as.character(prelimlag$dshortexp))
prelimlag$reserves <- as.numeric(as.character(prelimlag$reserves))
prelimlag$curraccgdp <- as.numeric(as.character(prelimlag$curraccgdp))
prelimlag$fdigdp <- as.numeric(as.character(prelimlag$fdigdp))
prelimlag$inflation <- as.numeric(as.character(prelimlag$inflation))
prelimlag$unsc <- as.numeric(as.character(prelimlag$unsc))
prelimlag$loansize <- as.numeric(as.character(prelimlag$loansize))
prelimlag$dur <- as.numeric(prelimlag$dur)
prelimlagimp <- mice(prelimlag, method = "cart", seed = 500, m = 1,
                  pred=quickpred(prelimlag, exclude=c("X", "ID", "ctry", "date", "type", 
                                                      "loansize")))
prelimlagimp <- complete(prelimlagimp)
write.csv(prelimlagimp, "prelimlagimp.csv")

####################### Load Data #############
## Note that before loading these datasets (as well as the text datasets used for text analysis in the
## following sections that are written directly to file from this script) I manually added the column 
## name "rowID" in Excel to the first column for each created dataset. These datasets when written to
## file from R include a rowID column without a name. The final datasets in the replication package
## already have this column name added (it is necessary to prevent merge errors). If you choose to 
## recreate the datasets using the above code, you will need to add the column name yourself.
prelim <- read.csv("prelim.csv")
prelimimp <- read.csv("prelimimp.csv")
prelimlag <- read.csv("prelimlag.csv")
prelimimplag <- read.csv("prelimlagimp.csv")
memsimplag <- read.csv("memsimp.csv")
####################### Data Pre-processing ##########

prelimlag$dur <- scale(as.numeric(as.character(prelimlag$dur))) # standardize all data to ease convergence of models
prelimlag$loansize <- log1p(as.numeric(as.character(prelimlag$loansize)))
prelimlag$quota <- scale(log1p(prelimlag$quota))
prelimlag$polity2 <- scale(prelimlag$polity2)
prelimlag$absidealimportantdiff <- scale(prelimlag$absidealimportantdiff)
prelimlag$absidealdiff <- scale(prelimlag$absidealdiff)
prelimlag$USaid <- scale(log1p(prelimlag$USaid))
prelimlag$DACaid <- scale(log1p(prelimlag$DACaid))
prelimlag$gdppc <- scale(prelimlag$gdppc)
prelimlag$gdp <- scale(log(prelimlag$gdp))
prelimlag$trade <- scale(prelimlag$trade)
prelimlag$reserves <- scale(log1p(prelimlag$reserves))
prelimlag$fdigdp <- scale(prelimlag$fdigdp)
prelimlag$dservexp <- scale(prelimlag$dservexp)
prelimlag$dshortexp <- scale(prelimlag$dshortexp)
prelimlag$ctry <- as.character(prelimlag$ctry)
prelimlag$ccode <- as.character(prelimlag$ccode)
prelimlag$ooamt <- scale(log1p(prelimlag$ooamt))
prelimlag$coopamt <- scale(log1p(prelimlag$coopamt))
prelimlag$inflation <- scale(prelimlag$inflation)
prelimlag$curraccgdp <- scale(prelimlag$curraccgdp)
prelimlag$prop_neoliberal <- scale(prelimlag$prop_neoliberal)
prelimlag$USaid[prelimlag$ctry %in% c("Portugal", "Ireland", "Cyprus")] <- 0 # code EU countries with missing data as no US aid
prelimlag$stone <- rowSums(prelimlag[,c("DEB","FIN","FP","EXT","RTP","SOE","LAB","PRI","SP","POV",
                                    "INS","ENV","OTH")]>0) # create number of categories variable

prelimimp$dur <- scale(as.numeric(as.character(prelimimp$dur))) # standardize all data to ease convergence of models
prelimimp$loansize <- log1p(as.numeric(as.character(prelimimp$loansize)))
prelimimp$quota <- scale(log1p(prelimimp$quota))
prelimimp$polity2 <- scale(prelimimp$polity2)
prelimimp$absidealimportantdiff <- scale(prelimimp$absidealimportantdiff)
prelimimp$absidealdiff <- scale(prelimimp$absidealdiff)
prelimimp$USaid <- scale(log1p(prelimimp$USaid))
prelimimp$DACaid <- scale(log1p(prelimimp$DACaid))
prelimimp$gdppc <- scale(prelimimp$gdppc)
prelimimp$gdp <- scale(log(prelimimp$gdp))
prelimimp$trade <- scale(prelimimp$trade)
prelimimp$reserves <- scale(log1p(prelimimp$reserves))
prelimimp$fdigdp <- scale(prelimimp$fdigdp)
prelimimp$dservexp <- scale(prelimimp$dservexp)
prelimimp$dshortexp <- scale(prelimimp$dshortexp)
prelimimp$ctry <- as.character(prelimimp$ctry)
prelimimp$ccode <- as.character(prelimimp$ccode)
prelimimp$ooamt <- scale(log1p(prelimimp$ooamt))
prelimimp$coopamt <- scale(log1p(prelimimp$coopamt))
prelimimp$inflation <- scale(prelimimp$inflation)
prelimimp$curraccgdp <- scale(prelimimp$curraccgdp)
prelimimp$prop_neoliberal <- scale(prelimimp$prop_neoliberal)
prelimimp$stone <- rowSums(prelimimp[,c("DEB","FIN","FP","EXT","RTP","SOE","LAB","PRI","SP","POV",
                                        "INS","ENV","OTH")]>0) # create number of categories variable

prelimimplag$dur <- scale(as.numeric(as.character(prelimimplag$dur))) # standardize data
prelimimplag$loansize <- log1p(as.numeric(as.character(prelimimplag$loansize)))
prelimimplag$quota <- scale(log1p(prelimimplag$quota))
prelimimplag$polity2 <- scale(prelimimplag$polity2)
prelimimplag$absidealimportantdiff <- scale(prelimimplag$absidealimportantdiff)
prelimimplag$absidealdiff <- scale(prelimimplag$absidealdiff)
prelimimplag$USaid <- scale(log1p(prelimimplag$USaid))
prelimimplag$DACaid <- scale(log1p(prelimimplag$DACaid))
prelimimplag$gdppc <- scale(prelimimplag$gdppc)
prelimimplag$gdp <- scale(log(prelimimplag$gdp))
prelimimplag$trade <- scale(prelimimplag$trade)
prelimimplag$reserves <- scale(log1p(prelimimplag$reserves))
prelimimplag$fdigdp <- scale(prelimimplag$fdigdp)
prelimimplag$dservexp <- scale(prelimimplag$dservexp)
prelimimplag$dshortexp <- scale(prelimimplag$dshortexp)
prelimimplag$ctry <- as.character(prelimimplag$ctry)
prelimimplag$ccode <- as.character(prelimimplag$ccode)
prelimimplag$ooamt <- scale(log1p(prelimimplag$ooamt))
prelimimplag$coopamt <- scale(log1p(prelimimplag$coopamt))
prelimimplag$inflation <- scale(prelimimplag$inflation)
prelimimplag$curraccgdp <- scale(prelimimplag$curraccgdp)
prelimimplag$prop_neoliberal <- scale(prelimimplag$prop_neoliberal)
prelimimplag$stone <- rowSums(prelimimplag[,c("DEB","FIN","FP","EXT","RTP","SOE","LAB","PRI","SP","POV",
                                        "INS","ENV","OTH")]>0) # create number of categories variable
prelimimplag2 <- prelimimplag[prelimimplag$PRGF == 0,] # robustness dataset excluding PRGF

memsimplag$quota <- scale(log1p(memsimplag$quota)) # standardize data
memsimplag$polity2 <- scale(memsimplag$polity2)
memsimplag$absidealimportantdiff <- scale(memsimplag$absidealimportantdiff)
memsimplag$absidealdiff <- scale(memsimplag$absidealdiff)
memsimplag$USaid <- scale(log1p(memsimplag$USaid))
memsimplag$DACaid <- scale(log1p(memsimplag$DACaid))
memsimplag$gdppc <- scale(memsimplag$gdppc)
memsimplag$gdp <- scale(log1p(memsimplag$gdp))
memsimplag$trade <- scale(memsimplag$trade)
memsimplag$reserves <- scale(log1p(memsimplag$reserves))
memsimplag$fdigdp <- scale(memsimplag$fdigdp)
memsimplag$dservexp <- scale(memsimplag$dservexp)
memsimplag$dshortexp <- scale(memsimplag$dshortexp)
memsimplag$ctry <- as.character(memsimplag$ctry)
memsimplag$ccode <- as.character(memsimplag$ccode)
memsimplag$inflation <- scale(memsimplag$inflation)
memsimplag$curraccgdp <- scale(memsimplag$curraccgdp)
memsimplag$coopamt <- scale(log1p(memsimplag$coopamt))
memsimplag$ooamt <- scale(log1p(memsimplag$ooamt))
tempimf <- aggregate(memsimplag$part, by = list(memsimplag$year), FUN = mean) # create dataset with ratio of member countries receiving support in a year
# this serves as a measure of budget constraint for the IMF
colnames(tempimf) <- c("year", "ratio") # change column names
tempc <- aggregate(memsimplag$part, by = list(memsimplag$ccode), FUN = mean) # create dataset with ratio of years that each country received support for IMF
# this serves as a measure of the likelihood a country will receive support in a year
colnames(tempc) <- c("ccode", "ratio") # change column names
memsimplag$imfratio <- NA # create excludable interaction variable for instrumentation
memsimplag$cratio <- NA # first create empty columns
for (i in 1:6398) {
  for (j in 1:37) {
    memsimplag$imfratio[i] <- ifelse(memsimplag$year[i] == tempimf$year[j], tempimf$ratio[j], 
                                     memsimplag$imfratio[i]) 
  }
} # then write ratios to the selection dataset
for (i in 1:6398) {
  for (j in 1:187) {
    memsimplag$cratio[i] <- ifelse(memsimplag$ccode[i] == tempc$ccode[j], tempc$ratio[j], 
                                     memsimplag$cratio[i]) 
  }
}
memsimplag$cratio <- scale(memsimplag$cratio) # standardize the ratios to ease convergence
memsimplag$imfratio <- scale(memsimplag$imfratio)
memsimplag$inter <- memsimplag$cratio*memsimplag$imfratio # interact the two ratios

####################### Main Tests ##################################################

ctryfe <- list(c("Country fixed effects", "Yes", "Yes")) # create rows denoting model specifications for regression tables
ctryre <- list(c("Country random effects", "Yes", "Yes"))
yearfe <- list(c("Year fixed effects", "Yes", "Yes"))
regionfe <- list(c("Region fixed effects", "Yes", "Yes"))
ipw <- list(c("Inverse probabiliy weights", "Yes", "Yes"))
nb <- list(c("Model type", "Negative binomial", "Negative binomial"))
pos <- list(c("Model type", "Poisson", "Poisson"))
prgf <- list(c("PRGF random effects", "Yes", "Yes"))

covs_full <- c("Outside option member", "Cooperative member", "Outside option amount", "Cooperative amount", 
               "Duration", "Quota", "PRGF", "Time from last IMF program", "Neoliberal ideology", 
               "Polity2", "UN voting (ideal pt dist from U.S.)", "UNSC member", "Checks", 
               "U.S. aid", "FDI / GDP", "GDPPC", "Openness", "Debt service / exports", 
               "Short-term debt / exports", "War", "Election year", "Year") # covariate labels

covs_bin <- c("Outside option member", "Cooperative member", "Outside option amount", 
              "Cooperative amount", "Year") # covariate labels for parsimonious models

# No controls models with lag and no imputation

regbiv <- glm.nb(count ~ oomem + coopmem + ooamt + coopamt + 
                 year + factor(ccode), data = prelimlag, control=glm.control(maxit=20), 
               na.action = "na.exclude")
(seregbiv <- coeftest(regbiv, vcov = vcovHC(regbiv, type = "HC0", cluster = prelimlag$ccode))) # parsimonious model count DV

regbiva <- glm.nb(stone ~ oomem + coopmem + ooamt + coopamt +
                 year + factor(ccode), data = prelimlag, control=glm.control(maxit=20))
(seregbiva <- coeftest(regbiva, vcov = vcovHC(regbiva, type = "HC0", cluster = prelimlag$ccode))) # parsimonious model categories DV

stargazer(regbiv, regbiva, se = list(seregbiv[,2], seregbiva[,2]), 
          dep.var.labels = c("Number of conditions", "Number of categories"),
          covariate.labels = covs_bin,
          style = "ajps", omit = c("ccode", "Constant"), add.lines = c(ctryfe, nb),
          omit.stat = c("ll", "aic", "theta")) # Table 1 in paper

# Main models with one year lags, imputation, ctry FEs, full controls and time trend
reg1 <- glm.nb(count ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                    unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec +
                    year + factor(ccode), data = prelimimplag, control=glm.control(maxit=20), 
               na.action = "na.exclude")
(sereg1 <- coeftest(reg1, vcov = vcovHC(reg1, type = "HC0", cluster = prelimimplag$ccode))) # Full nbreg model count DV 

reg2 <- glm.nb(stone ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                 unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + 
                 year + factor(ccode), data = prelimimplag, control=glm.control(maxit=20))
(sereg2 <- coeftest(reg2, vcov = vcovHC(reg2, type = "HC0", cluster = prelimimplag$ccode))) # Full nbreg model categories DV

stargazer(reg1, reg2, se = list(sereg1[,2], sereg2[,2]), 
          dep.var.labels = c("Number of conditions", "Number of categories"),
          covariate.labels = covs_full,
          style = "ajps", omit = c("ccode", "Constant"), add.lines = c(ctryfe, nb),
          omit.stat = c("ll", "aic", "theta")) # Table 2 in paper

####################### Selection Adjustment ##########
c2 <- c("Budget constraint:Participation rate", "Outside option member", "Cooperative member", 
        "Outside option amount", "Cooperative amount", "Quota", "Time to last IMF program", "Polity2", 
        "Reserves", "GDPPC", "Current account / GDP", "UNSC member", "U.S. aid", "DAC aid", 
        "UN voting", "FDI / GDP", "Inflation", "Openness", "Debt service / exports", 
        "Short-term debt / exports", "War", "Election year") # covariate labels for first stage
prob <- list(c("Model type", "Probit"))
## First stage selection model

partreg <- glm(part ~ inter + oomem + coopmem + ooamt + coopamt + quota + last + polity2 + reserves + gdppc + curraccgdp +
                 unsc + USaid + DACaid + absidealdiff + fdigdp + inflation + trade + dservexp + dshortexp +
                 war + elec + factor(year) + factor(ccode), family = quasibinomial(link = "probit"), 
               data = memsimplag, control=glm.control(maxit=200))
(separt <- coeftest(partreg, vcov = vcovHC(partreg, type = "HC0", cluster = memsimplag$ccode))) # binomial probit with instrument

partregsub <- glm(part ~ oomem + coopmem + ooamt + coopamt + quota + last + polity2 + reserves + gdppc + curraccgdp +
                 unsc + USaid + DACaid + absidealdiff + fdigdp + inflation + trade + dservexp + dshortexp +
                 war + elec + factor(year) + factor(ccode), family = quasibinomial(link = "probit"), 
               data = memsimplag, control=glm.control(maxit=200))
(separt2 <- coeftest(partregsub, vcov = vcovHC(partregsub, type = "HC0", cluster = memsimplag$ccode))) # delete instrument for weak instrument test

memsimplag <- memsimplag[!is.na(memsimplag$inter),] # eliminate observations with NA interaction
memsimplag$fit <- fitted(partreg) # add predicted probabilities to data frame

## Merge predicted probabilities into selection data

fitdat <- cbind(as.character(memsimplag$ccode), memsimplag$year, memsimplag$fit) # extract probabilities
fitdat <- unique(fitdat) # eliminate duplicates
colnames(fitdat) <- c("ccode", "year", "fit") # rename columns
prelimimplag4 <- merge(prelimimplag, fitdat, c("ccode", "year")) # merge data
prelimimplag4 <- unique(prelimimplag4) # eliminate duplicates
prelimimplag4$fit <- as.numeric(as.character(prelimimplag4$fit)) # change data format

stargazer(partreg, se = list(separt[,2]),
          dep.var.labels = c("Participation in IMF program"),
          covariate.labels = c2,
          style = "ajps", omit = c("ccode", "year", "Constant"), 
          add.lines = c(ctryfe, yearfe, prob),
          omit.stat = c("ll", "aic", "theta")) # Table A13 in appendix

## Second stage with inverse probability weights

regfa <- glm(count ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + 
               polity2 + absidealdiff + unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + 
               dshortexp + war + elec + year + factor(ctry), 
             data = prelimimplag4, weights = fit, family = poisson(link = "log"), control=glm.control(maxit=200))
(seregfa <- coeftest(regfa, vcov = vcovHC(regfa, type = "HC0", cluster = prelimimplag4$ctry))) # Count DV with IPW

regfb <- glm(stone ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + 
               polity2 + absidealdiff + unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + 
               dshortexp + war + elec + year + factor(ctry), 
             data = prelimimplag4, weights = fit, family = poisson(link = "log"), control=glm.control(maxit=200))
(seregfb <- coeftest(regfb, vcov = vcovHC(regfb, type = "HC0", cluster = prelimimplag4$ctry))) # Categories DV with IPW

stargazer(regfa, regfb, se = list(seregfa[,2], seregfb[,2]),
          dep.var.labels = c("Number of prior actions", "Number of categories"),
          covariate.labels = covs_full,
          style = "ajps", omit = c("ctry","Constant"),
          add.lines = c(ctryfe, ipw, pos),
          omit.stat = c("ll", "aic", "theta")) # Table A14 in appendix

## Weak instrument test
anova(partreg, partregsub, test = "F") # F of 30 tells us the difference is significant

## Heckman selection

selsub <- cbind(as.character(prelimimplag$ccode), prelimimplag$year, prelimimplag$count,  
                prelimimplag$stone, prelimimplag$dur, prelimimplag$PRGF, prelimimplag$checks,
                prelimimplag$prop_neoliberal) # add variables not currently in the selection data to the selection data
colnames(selsub) <- c("ccode","year","count", "stone", "dur", "PRGF", "checks", "prop_neoliberal") # change column names
selsub <- data.frame(selsub) # convert to frame
selsub$count <- as.numeric(as.character(selsub$count)) # change data format
selsub$stone <- as.numeric(as.character(selsub$stone))
selsub$dur <- as.numeric(as.character(selsub$dur))
selsub$PRGF <- as.numeric(as.character(selsub$PRGF))
selsub$checks <- as.numeric(as.character(selsub$checks))
selsub$prop_neoliberal <- as.numeric(as.character(selsub$prop_neoliberal))
selsub <- aggregate(list(selsub$count, selsub$stone, selsub$dur, selsub$PRGF, selsub$checks, selsub$prop_neoliberal), 
                    by = list(selsub$ccode, selsub$year), FUN = max) # aggregate data by country-year
colnames(selsub) <- c("ccode","year","count", "stone", "dur", "PRGF", "checks", "prop_neoliberal") # rename columns
sel <- merge(unique(memsimplag), unique(selsub), by = c("ccode", "year"), all.x = T) # merge data
sel <- unique(sel) # eliminate duplicates
sel$part <- as.logical(sel$part) # change to logical for Heckman
sel2 <- sel[sel$part == 1,] # subset to be used in second stage

regfc <- selection(selection = part ~ inter + oomem + coopmem + ooamt + coopamt + quota + last + polity2 + gdppc +
                     unsc + USaid + absidealdiff + fdigdp + inflation + trade + dservexp + dshortexp +
                     war + elec + factor(year) + factor(ccode), 
                outcome = count ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + unsc + checks +
                     USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + year + factor(ccode), 
                data = sel) # Heckman model
summary(regfc) # will not converge, as noted in on appendix p. 29

regfd <- heckit(selection = part ~ inter + oomem + coopmem + ooamt + coopamt + quota + last + polity2 + gdppc +
                  unsc + USaid + absidealdiff + fdigdp + inflation + trade + dservexp + dshortexp +
                  war + elec + factor(ccode) + factor(year), 
                outcome = stone ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + unsc + checks +
                  USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + year + factor(ccode), 
                data = sel, inst = ~ inter) # Heckman model
summary(regfd) # will not converge, as noted in on appendix p. 29

####################### Other Robustness Checks ################
covs_full_noimp <- c("Outside option member", "Cooperative member", "Outside option amount", "Cooperative amount", 
               "Duration", "Quota", "PRGF", "Time from last IMF program", 
               "Polity2", "UN voting (ideal pt dist from U.S.)", "UNSC member", "Checks", 
               "U.S. aid", "FDI / GDP", "GDPPC", "Openness", "War", "Election year", "Year") # covariate labels excluding debt service, short-term debt, and prop neoliberal to conserve observations
## No imputation all covariates (ctry FEs, full controls)
reg4 <- glm.nb(count ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + polity2 + absidealdiff + 
                    unsc + checks + USaid + fdigdp + gdppc + trade + war + elec + 
                    year + factor(ctry), data = prelimlag, control=glm.control(maxit=20))
(sereg4 <- coeftest(reg4, vcov = vcovHC(reg4, type = "HC0", cluster = prelimlag$ctry))) # Count DV

reg4a <- glm.nb(stone ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + polity2 + absidealdiff + unsc + 
                  checks + USaid + fdigdp + gdppc + trade + war + elec + 
                  year + factor(ctry), data = prelimlag, control=glm.control(maxit=20))
(sereg4a <- coeftest(reg4a, vcov = vcovHC(reg4a, type = "HC0", cluster = prelimlag$ctry))) # Categories DV

stargazer(reg4,reg4a, se = list(sereg4[,2], sereg4a[,2]),
          dep.var.labels = c("Number of conditions", "Number of categories"),
          covariate.labels = covs_full_noimp,
          style = "ajps", omit = c("ctry", "Constant"), add.lines = c(ctryfe, nb),
          omit.stat = c("ll", "aic", "theta")) # Table A16 in appendix

## No lag (imputation, ctry FEs, full controls)

reg5 <- glm.nb(count ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                 unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + 
                 year + factor(ctry), data = prelimimp, control=glm.control(maxit=20))
(sereg5 <- coeftest(reg5, vcov = vcovHC(reg5, type = "HC0", cluster = prelimimp$ctry))) # count DV

reg5a <- glm.nb(stone ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + unsc + 
               checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + 
               year + factor(ctry), data = prelimimp, family = poisson(link = "log"), control=glm.control(maxit=20))
(sereg5a <- coeftest(reg5a, vcov = vcovHC(reg5a, type = "HC0", cluster = prelimimp$ctry))) # categories DV

stargazer(reg5,reg5a, se = list(sereg5[,2], sereg5a[,2]),
          dep.var.labels = c("Number of conditions", "Number of categories"),
          covariate.labels = covs_full,
          style = "ajps", omit = c("ctry"), add.lines = c(ctryfe, nb),
          omit.stat = c("ll", "aic", "theta")) # Table A17 in appendix

## No imputation (ctry FEs, limited controls)

covs_limit <- c("Outside option member", "Cooperative member", "Outside option amount", "Cooperative amount", 
               "UN voting (ideal pt dist from U.S.)", "UNSC member", 
               "U.S. aid", "Year") # modified covariate labels

reg6 <- glm.nb(count ~ oomem + coopmem + ooamt + coopamt + absidealdiff + 
                 unsc + USaid + year + factor(ctry), data = prelimlag, control=glm.control(maxit=20))
(sereg6 <- coeftest(reg6, vcov = vcovHC(reg6, type = "HC0", cluster = prelimlag$ctry))) # count DV

reg6a <- glm.nb(stone ~ oomem + coopmem + ooamt + coopamt + absidealdiff + 
               unsc + USaid + year + factor(ctry), data = prelimlag, control=glm.control(maxit=20))
(sereg6a <- coeftest(reg6a, vcov = vcovHC(reg6a, type = "HC0", cluster = prelimlag$ctry))) # categories DV

stargazer(reg6,reg6a, se = list(sereg6[,2], sereg6a[,2]),
          dep.var.labels = c("Number of conditions", "Number of categories"),
          covariate.labels = covs_limit,
          style = "ajps", omit = c("ctry", "Constant"), add.lines = c(ctryfe, nb),
          omit.stat = c("ll", "aic", "theta")) # Table A15 in appendix

## Drop Greece 

prelimimplagg <- prelimimplag[!prelimimplag$ctry == "Greece",] # exclude Greece from sample

reg7 <- glm.nb(count ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                 unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec +
                 year + factor(ctry), data = prelimimplagg, control=glm.control(maxit=20), 
               na.action = "na.exclude") # count DV
(sereg7 <- coeftest(reg7, vcov = vcovHC(reg7, type = "HC0", cluster = prelimimplagg$ctry))) 

reg7a <- glm.nb(stone ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                 unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec +
                 year + factor(ctry), data = prelimimplagg, control=glm.control(maxit=20), 
               na.action = "na.exclude") # categories DV
(sereg7a <- coeftest(reg7a, vcov = vcovHC(reg7a, type = "HC0", cluster = prelimimplagg$ctry)))

stargazer(reg7, reg7a, se = list(sereg7[,2], sereg7a[,2]), 
          dep.var.labels = c("Number of conditions", "Number of categories"),
          covariate.labels = covs_full,
          style = "ajps", omit = c("ctry", "Constant"), add.lines = c(ctryfe, nb),
          omit.stat = c("ll", "aic", "theta")) # Table A19 in appendix

## Year FEs instead of time trend

reg9 <- glm.nb(count ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                 unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec +
                 factor(ccode) + factor(year), data = prelimimplag, control=glm.control(maxit=20), 
               na.action = "na.exclude")
(sereg9 <- coeftest(reg9, vcov = vcovHC(reg9, type = "HC0", cluster = prelimimplag$ccode))) # count DV 

reg9a <- glm.nb(stone ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                  unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + 
                  factor(ccode) + factor(year), data = prelimimplag, control=glm.control(maxit=20))
(sereg9a <- coeftest(reg9a, vcov = vcovHC(reg9a, type = "HC0", cluster = prelimimplag$ccode))) # categories DV

stargazer(reg9, reg9a, se = list(sereg9[,2], sereg9a[,2]), 
          dep.var.labels = c("Number of conditions", "Number of categories"),
          covariate.labels = covs_full,
          style = "ajps", omit = c("ccode", "year", "Constant"), add.lines = c(ctryfe, yearfe, nb),
          omit.stat = c("ll", "aic", "theta")) # Table A21 in appendix

## Info sharing and Co-financing separately

covs_r3 <- c("Outside option member", "Co-financing", "Information Sharing", "Outside option amount", "Cooperative amount", 
             "Duration", "Quota", "PRGF", "Time from last IMF program", "Neoliberal ideology", 
             "Polity2", "UN voting (ideal pt dist from U.S.)", "UNSC member", "Checks", 
             "U.S. aid", "FDI / GDP", "GDPPC", "Openness", "Debt service / exports", 
             "Short-term debt / exports", "War", "Election year", "Year") # adjusted covariate labels

reg11 <- glm.nb(count ~ oomem + cofin + infoshare + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                  unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + year +
                  factor(ccode), data = prelimimplag, control=glm.control(maxit=20), 
                na.action = "na.exclude")
(sereg11 <- coeftest(reg11, vcov = vcovHC(reg11, type = "HC0", cluster = prelimimplag$ccode))) # count DV 

reg11a <- glm.nb(stone ~ oomem + cofin + infoshare + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                  unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + year +
                  factor(ccode), data = prelimimplag, control=glm.control(maxit=20), 
                na.action = "na.exclude")
(sereg11a <- coeftest(reg11a, vcov = vcovHC(reg11a, type = "HC0", cluster = prelimimplag$ccode))) # categories DV 

stargazer(reg11, reg11a, se = list(sereg11[,2], sereg11a[,2]), 
          dep.var.labels = c("Number of conditions", "Number of categories"),
          covariate.labels = covs_r3,
          style = "ajps", omit = c("ccode", "Constant"), add.lines = c(ctryfe, nb),
          omit.stat = c("ll", "aic", "theta")) # Table A24 in appendix

## Waiver adjustment robustness check

prelimimplag$countw <- prelimimplag$count - prelimimplag$waiver # alternate count measure where waived conditions are subtracted from total
reg12 <- glm.nb(countw ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                  unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + year +
                  factor(ccode), data = prelimimplag, control=glm.control(maxit=20), 
                na.action = "na.exclude")
(sereg12 <- coeftest(reg12, vcov = vcovHC(reg12, type = "HC0", cluster = prelimimplag$ccode))) # adjusted count DV 

stargazer(reg12, se = list(sereg12[,2]), 
          dep.var.labels = c("Number of conditions"),
          covariate.labels = covs_full,
          style = "ajps", omit = c("ccode", "Constant"), add.lines = c(ctryfe, nb),
          omit.stat = c("ll", "aic", "theta")) # Table A27

## Drop ESM

esmless <- prelimimplag[prelimimplag$ESMm != 1,] # exclude ESM members from analysis

reg15 <- glm.nb(count ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                 unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec +
                 year + factor(ccode), data = esmless, control=glm.control(maxit=20), 
               na.action = "na.exclude")
(sereg15 <- coeftest(reg15, vcov = vcovHC(reg15, type = "HC0", cluster = esmless$ccode))) # Full nbreg model count DV 

reg16 <- glm.nb(stone ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                 unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + 
                 year + factor(ccode), data = esmless, control=glm.control(maxit=20))
(sereg16 <- coeftest(reg16, vcov = vcovHC(reg16, type = "HC0", cluster = esmless$ccode))) # Full nbreg model categories DV

stargazer(reg15, reg16, se = list(sereg15[,2], sereg16[,2]), 
          dep.var.labels = c("Number of conditions", "Number of categories"),
          covariate.labels = covs_full,
          style = "ajps", omit = c("ccode", "Constant"), add.lines = c(ctryfe, nb),
          omit.stat = c("ll", "aic", "theta")) # Table A20 in appendix

## Deference control (by controlling for CMIM membership explicitly)

covs_r4 <- c("Outside option member", "Cooperative member", "Outside option amount", "Cooperative amount", 
             "Duration", "Quota", "PRGF", "Time from last IMF program", "Neoliberal ideology", 
             "Polity2", "UN voting (ideal pt dist from U.S.)", "UNSC member", "Checks", 
             "U.S. aid", "FDI / GDP", "GDPPC", "Openness", "Debt service / exports", 
             "Short-term debt / exports", "War", "Election year", "Deference", "Year") # adjusted covariate labels

reg17 <- glm.nb(count ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                  unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec +
                  year + CMIMm + factor(ccode), data = prelimimplag, control=glm.control(maxit=20), 
                na.action = "na.exclude")
(sereg17 <- coeftest(reg17, vcov = vcovHC(reg17, type = "HC0", cluster = prelimimplag$ccode))) # Full nbreg model count DV 

reg18 <- glm.nb(stone ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                  unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + 
                  year + CMIMm + factor(ccode), data = prelimimplag, control=glm.control(maxit=20))
(sereg18 <- coeftest(reg18, vcov = vcovHC(reg18, type = "HC0", cluster = prelimimplag$ccode))) # Full nbreg model categories DV

stargazer(reg17, reg18, se = list(sereg17[,2], sereg18[,2]), 
          dep.var.labels = c("Number of conditions", "Number of categories"),
          covariate.labels = covs_r4,
          style = "ajps", omit = c("ccode"), add.lines = c(ctryfe, nb),
          omit.stat = c("ll", "aic", "theta")) # Table A22 in appendix

## Bilateral bailouts and currency swaps

stdata <- read.table("schneider_tobin_bailouts_final.tab", fill=T) # load bilateral bailout data from Schenider and Tobin (2020) IO https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/HJ9SWF
header.true <- function(df) {
  names(df) <- as.character(unlist(df[1,]))
  df[-1,]
} # function to make first row the header names
stdata <- header.true(stdata)
stdata[stdata == ""] <- NA # write blanks as NA
stdata <- aggregate(as.numeric(stdata$bilateral_bailout), by = list(stdata$crisis_country, stdata$year), 
                    FUN = sum) # aggregate bilateral bailouts by recipient country-year
colnames(stdata) <- c("ctry", "year", "bilateral") # rename cols
stdata$ccode <- countrycode(stdata$ctry, origin = "country.name", destination = "iso3c") # create country code variable
stdata$year <- as.numeric(stdata$year) + 1 # lag year
stcheck <- merge(unique(prelimimplag), unique(stdata), by = c("ccode", "year"), all.x=T) # merge with main data
stcheck$bilateral[is.na(stcheck$bilateral)] <- 0 # recode NA bilateral as 0
stcheck$bilateral <- scale(as.numeric(stcheck$bilateral)) # scale bilateral to ease convergence
stcheck <- stcheck[stcheck$year < 2011,] # restrict to when bilateral data available

bsa <- read.csv("bsa.csv") # load hand-coded data on bilateral swap arragements with US, China, and Japan for IMF participant countries
bsa$ccode <- countrycode(bsa$ctry, origin = "country.name", destination = "iso3c") # create ccodes
stcheck <- merge(unique(stcheck), unique(bsa), by = c("ccode", "year"), all.x=T) # merge
stcheck$numswap[is.na(stcheck$numswap)] <- 0 # mark NA swaps as 0

covs_r5 <- c("Outside option member", "Cooperative member", "Outside option amount", "Cooperative amount", 
             "Duration", "Quota", "PRGF", "Time from last IMF program", "Neoliberal ideology", 
             "Polity2", "UN voting (ideal pt dist from U.S.)", "UNSC member", "Checks", 
             "U.S. aid", "FDI / GDP", "GDPPC", "Openness", "Debt service / exports", 
             "Short-term debt / exports", "War", "Election year", "Bilateral rescues", "Bilateral swaps", "Year") # adjusted covariate labels

reg19 <- glm.nb(count ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                  unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec +
                  year + bilateral + numswap + factor(ccode), data = stcheck, control=glm.control(maxit=20), 
                na.action = "na.exclude")
(sereg19 <- coeftest(reg19, vcov = vcovHC(reg19, type = "HC0", cluster = stcheck$ccode))) # Full nbreg model count DV 

reg20 <- glm.nb(stone ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                  unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + 
                  year + bilateral + numswap + factor(ccode), data = stcheck, control=glm.control(maxit=20))
(sereg20 <- coeftest(reg20, vcov = vcovHC(reg20, type = "HC0", cluster = stcheck$ccode))) # Full nbreg model categories DV

stargazer(reg19, reg20, se = list(sereg19[,2], sereg20[,2]), 
          dep.var.labels = c("Number of conditions", "Number of categories"),
          covariate.labels = covs_r5,
          style = "ajps", omit = "ccode", add.lines = c(ctryfe, nb),
          omit.stat = c("ll", "aic", "theta")) # Table A23 in appendix

## Weak states for exogenous cooperation ##

quantile(prelimimplag$gdp)
weakstates <- data.frame(prelimimplag[prelimimplag$gdp < 0.55,])

reg21 <- glm.nb(count ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                 unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec +
                 year + factor(ccode), data = weakstates, control=glm.control(maxit=20), 
               na.action = "na.exclude")
(sereg21 <- coeftest(reg21, vcov = vcovHC(reg21, type = "HC0", cluster = weakstates$ccode))) # Full nbreg model count DV 

reg22 <- glm.nb(stone ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                 unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + 
                 year + factor(ccode), data = weakstates, control=glm.control(maxit=20))
(sereg22 <- coeftest(reg22, vcov = vcovHC(reg22, type = "HC0", cluster = weakstates$ccode))) # Full nbreg model categories DV

stargazer(reg21, reg22, se = list(sereg21[,2], sereg22[,2]), 
          dep.var.labels = c("Number of conditions", "Number of categories"),
          covariate.labels = covs_full,
          style = "ajps", omit = "ccode", add.lines = c(ctryfe, nb),
          omit.stat = c("ll", "aic", "theta")) # Table A12 in appendix

## WB loans control
wb <- read.csv("wbloans.csv") # load World Bank loan data (IBRD loans and IDA credits) from WDI downloaded August 2021
wb <- wb[-c(2,3)] # delete extra cols
colnames(wb) <- c("year", "ccode", "wbloans") # rename cols
wb$year <- as.numeric(wb$year) # format
wb$wbloans <- as.numeric(wb$wbloans) # format
wb$wbloans[is.na(wb$wbloans)] <- 0 # change NA to 0 because no loans
wbcheck <- merge(unique(prelimimplag), unique(wb), by = c("ccode", "year"), all.x=T) # merge with main data
wbcheck$wbloans <- scale(log1p(wbcheck$wbloans)) # standardize and log

reg23 <- glm.nb(count ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                  unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + wbloans +
                  year + factor(ccode), data = wbcheck, control=glm.control(maxit=20), 
                na.action = "na.exclude")
(sereg23 <- coeftest(reg23, vcov = vcovHC(reg23, type = "HC0", cluster = wbcheck$ccode))) # Full nbreg model count DV 

reg24 <- glm.nb(stone ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                  unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + wbloans +
                  year + factor(ccode), data = wbcheck, control=glm.control(maxit=20))
(sereg24 <- coeftest(reg24, vcov = vcovHC(reg24, type = "HC0", cluster = wbcheck$ccode))) # Full nbreg model categories DV

covs_r6 <- c("Outside option member", "Cooperative member", "Outside option amount", "Cooperative amount", 
             "Duration", "Quota", "PRGF", "Time from last IMF program", "Neoliberal ideology", 
             "Polity2", "UN voting (ideal pt dist from U.S.)", "UNSC member", "Checks", 
             "U.S. aid", "FDI / GDP", "GDPPC", "Openness", "Debt service / exports", 
             "Short-term debt / exports", "War", "Election year", "World Bank Loans", "Year") # adjusted covariate labels

stargazer(reg23, reg24, se = list(sereg23[,2], sereg24[,2]), 
          dep.var.labels = c("Number of conditions", "Number of categories"),
          covariate.labels = covs_r6,
          style = "ajps", omit = "ccode", add.lines = c(ctryfe, nb),
          omit.stat = c("ll", "aic", "theta")) # Table A26 in appendix

# Sensitivity analysis

reg1l <- lm(count ~ oomem + dur + quota + last + PRGF + prop_neoliberal + polity2 + absidealdiff + 
             unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec +
             year + factor(ccode), data = prelimimplag)
(sereg1l <- coeftest(reg1l, vcov = vcovHC(reg1l, type = "HC0", cluster = prelimimplag$ccode))) 
# I use lm here because sensemakr only takes linear regression inputs. Results are similar to the neg bin. I include
# my focal treatment variable (oomem) along with the full cohort of controls in line with the recommendations of the 
# package materials, and I test sensitivity for both DVs at the 0.1 significance level.

reg2l <- lm(stone ~ oomem + dur + quota + last + PRGF + prop_neoliberal + polity2 + absidealdiff + 
             unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec +
             year + factor(ccode), data = prelimimplag)
(sereg2l <- coeftest(reg2l, vcov = vcovHC(reg2l, type = "HC0", cluster = prelimimplag$ccode))) # lm with treatment for sensitivity 

sens <- sensemakr(model = reg1l, 
          treatment = "oomem",
          benchmark_covariates = "USaid",
          alpha = 0.1,
          kd = 1:3)

summary(sens)

ovb_minimal_reporting(sens, format = "latex") # Table A28 in appendix

sens <- sensemakr(model = reg2l, 
                  treatment = "oomem",
                  benchmark_covariates = "USaid",
                  alpha = 0.1,
                  kd = 1:3)

summary(sens)

ovb_minimal_reporting(sens, format = "latex") # Table A29 in appendix

# Multi-level modeling on PRGF
prelimimplag$PRGF <- as.factor(prelimimplag$PRGF)
prelimimplag$ccode <- as.factor(prelimimplag$ccode)
reg27 <- glmer.nb(count ~ oomem + coopmem + ooamt + coopamt + dur + quota + last + prop_neoliberal + polity2 + absidealdiff + 
                 unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec +
                 year + ccode + (1 | PRGF), data = prelimimplag,
                 control=glmerControl(optCtrl=list(maxfun=20)), nAGQ=0) # incorporates random intercept for PRGF
summary(reg27)

reg28 <- glmer.nb(stone ~ oomem + coopmem + ooamt + coopamt + dur + quota + last + prop_neoliberal + polity2 + absidealdiff + 
                 unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + 
                 year + ccode + (1 | PRGF), data = prelimimplag, 
                 control=glmerControl(optCtrl=list(maxfun=20)), nAGQ=0) # incorporates random intercept for PRGF
summary(reg28)

covs_r7 <- c("Outside option member", "Cooperative member", "Outside option amount", "Cooperative amount", 
             "Duration", "Quota", "Time from last IMF program", "Neoliberal ideology", 
             "Polity2", "UN voting (ideal pt dist from U.S.)", "UNSC member", "Checks", 
             "U.S. aid", "FDI / GDP", "GDPPC", "Openness", "Debt service / exports", 
             "Short-term debt / exports", "War", "Election year", "Year") # adjusted covariate labels

stargazer(reg27, reg28,
          dep.var.labels = c("Number of conditions", "Number of categories"),
          covariate.labels = covs_r7,
          style = "ajps", omit = "ccode", add.lines = c(ctryfe, prgf, nb),
          omit.stat = c("ll", "aic", "theta")) # Table A25 in appendix

# Region fixed effects

prelimimplag$region <- countrycode(prelimimplag$ccode, origin = "iso3c", destination = "region")
reg29 <- glm.nb(count ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                 unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec +
                 year + factor(ccode) + factor(region), data = prelimimplag, control=glm.control(maxit=20), 
               na.action = "na.exclude")
(sereg29 <- coeftest(reg29, vcov = vcovHC(reg29, type = "HC0", cluster = prelimimplag$ccode)))

reg30 <- glm.nb(stone ~ oomem + coopmem + ooamt + coopamt + dur + quota + PRGF + last + prop_neoliberal + polity2 + absidealdiff + 
                 unsc + checks + USaid + fdigdp + gdppc + trade + dservexp + dshortexp + war + elec + 
                 year + factor(ccode) + factor(region), data = prelimimplag, control=glm.control(maxit=20))
(sereg30 <- coeftest(reg30, vcov = vcovHC(reg30, type = "HC0", cluster = prelimimplag$ccode)))

stargazer(reg29, reg30,
          dep.var.labels = c("Number of conditions", "Number of categories"),
          covariate.labels = covs_full,
          style = "ajps", omit = c("ccode", "region"), add.lines = c(ctryfe, regionfe, nb),
          omit.stat = c("ll", "aic", "theta")) # Table A18 in appendix

####################### Robustness Plots ####
varnames <- c("Outside option member", "Cooperative member", 
              "Outside option member", "Cooperative member",
              "Outside option member", "Cooperative member", 
              "Outside option member", "Cooperative member",
              "Outside option member", "Cooperative member", 
              "Outside option member", 
              "Outside option member", "Cooperative member",
              "Outside option member", "Cooperative member", 
              "Outside option member", "Cooperative member", 
              "Outside option member", "Cooperative member",
              "Outside option member", "Cooperative member",
              "Outside option member", "Cooperative member", 
              "Outside option member", "Cooperative member",
              "Outside option member", "Cooperative member") # variable names to appear in plots
coef3 <- seregfa[,1]
coef3 <- unname(coef3[c(2,3)])
coef3a <- sereg5[,1]
coef3a <- unname(coef3a[c(2,3)])
coef3b <- sereg7[,1]
coef3b <- unname(coef3b[c(2,3)])
coef3c <- sereg9[,1]
coef3c <- unname(coef3c[c(2,3)])
coef3d <- sereg12[,1]
coef3d <- unname(coef3d[c(2,3)])
coef3e <- sereg11[,1]
coef3e <- unname(coef3e[c(2)])
coef3g <- sereg6[,1]
coef3g <- unname(coef3g[c(2,3)])
coef3h <- sereg15[,1]
coef3h <- unname(coef3h[c(2,3)])
coef3i <- sereg17[,1]
coef3i <- unname(coef3i[c(2,3)])
coef3j <- sereg19[,1]
coef3j <- unname(coef3j[c(2,3)])
coef3k <- sereg21[,1]
coef3k <- unname(coef3k[c(2,3)])
coef3l <- sereg23[,1]
coef3l <- unname(coef3l[c(2,3)])
coef3m <- summary(reg27)$coefficients
coef3m <- unname(coef3m[c(2,3)])
coef3n <- sereg29[,1]
coef3n <- unname(coef3n[c(2,3)])

se3 <- seregfa[,2]
se3 <- unname(se3[c(2,3)])
se3a <- sereg5[,2]
se3a <- unname(se3a[c(2,3)])
se3b <- sereg7[,2]
se3b <- unname(se3b[c(2,3)])
se3c <- sereg9[,2]
se3c <- unname(se3c[c(2,3)])
se3d <- sereg12[,2]
se3d <- unname(se3d[c(2,3)])
se3e <- sereg11[,2]
se3e <- unname(se3e[c(2)])
se3g <- sereg6[,2]
se3g <- unname(se3g[c(2,3)])
se3h <- sereg15[,2]
se3h <- unname(se3h[c(2,3)])
se3i <- sereg17[,2]
se3i <- unname(se3i[c(2,3)])
se3j <- sereg19[,2]
se3j <- unname(se3j[c(2,3)])
se3k <- sereg21[,2]
se3k <- unname(se3k[c(2,3)])
se3l <- sereg23[,2]
se3l <- unname(se3l[c(2,3)])
se3m <- sqrt(diag(vcov(reg27)))
se3m <- unname(se3m[c(2,3)])
se3n <- sereg29[,2]
se3n <- unname(se3n[c(2,3)])

low3 <- coef3 - (1.645*se3)
high3 <- coef3 + (1.645*se3)
low3a <- coef3a - (1.645*se3a)
high3a <- coef3a + (1.645*se3a)
low3b <- coef3b - (1.645*se3b)
high3b <- coef3b + (1.645*se3b)
low3c <- coef3c - (1.645*se3c)
high3c <- coef3c + (1.645*se3c)
low3d <- coef3d - (1.645*se3d)
high3d <- coef3d + (1.645*se3d)
low3e <- coef3e - (1.645*se3e)
high3e <- coef3e + (1.645*se3e)
low3g <- coef3g - (1.645*se3g)
high3g <- coef3g + (1.645*se3g)
low3h <- coef3h - (1.645*se3h)
high3h <- coef3h + (1.645*se3h)
low3i <- coef3i - (1.645*se3i)
high3i <- coef3i + (1.645*se3i)
low3j <- coef3j - (1.645*se3j)
high3j <- coef3j + (1.645*se3j)
low3k <- coef3k - (1.645*se3k)
high3k <- coef3k + (1.645*se3k)
low3l <- coef3l - (1.645*se3l)
high3l <- coef3l + (1.645*se3l)
low3m <- coef3m - (1.645*se3m)
high3m <- coef3m + (1.645*se3m)
low3n <- coef3n - (1.645*se3n)
high3n <- coef3n + (1.645*se3n)

resModel3 <- data.frame(param = varnames, 
                        low2.5 = c(low3, low3a, low3b, low3c, low3d, low3e, low3g, low3h, low3i, low3j,
                                   low3k, low3l, low3m, low3n), 
                        est = c(coef3, coef3a, coef3b, coef3c, coef3d, coef3e, coef3g, coef3h, coef3i, coef3j,
                                coef3k, coef3l, coef3m, coef3n), 
                        up2.5 = c(high3, high3a, high3b, high3c, high3d, high3e, high3g, high3h, high3i, high3j, 
                                  high3k, high3l, high3m, high3n)) 
resModel3$Model <- modelnames <- c("Selection", "Selection", 
                                   "No lag", "No lag", 
                                   "Drop Greece", "Drop Greece", 
                                   "Year FEs", "Year FEs", 
                                   "Drop waivers", "Drop waivers",
                                   "Disagg cooperation",  
                                   "No imputation", "No imputation",
                                   "Drop ESM", "Drop ESM",
                                   "Deference control", "Deference control",
                                   "Bilateral controls", "Bilateral controls",
                                   "Weak states", "Weak states",
                                   "World Bank control", "World Bank control",
                                   "Mutlilevel", "Multilevel",
                                   "Region FEs", "Region FEs")

limits <- aes(ymax  = up2.5, ymin = low2.5) 

plot3 <- ggplot(resModel3, aes(x=param, y=est, color = Model)) + 
  geom_pointrange(limits, position=position_dodge(width=0.2)) + 
  geom_hline(yintercept=0, color="red") + 
  theme_bw() + 
  theme(text = element_text(size=22), axis.title.y = element_blank(), 
        legend.position = "none") +
  ylab(NULL) + 
  coord_flip() +
  ggtitle("Conditions") +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_x_discrete(limits = rev(levels(as.factor(resModel3$param))))

varnames2 <- c("Outside option member", "Cooperative member", 
              "Outside option member", "Cooperative member", 
              "Outside option member", "Cooperative member", 
              "Outside option member", "Cooperative member", 
              "Outside option member",  
              "Outside option member", "Cooperative member", 
              "Outside option member", "Cooperative member",
              "Outside option member", "Cooperative member", 
              "Outside option member", "Cooperative member",
              "Outside option member", "Cooperative member",
              "Outside option member", "Cooperative member",
              "Outside option member", "Cooperative member", 
              "Outside option member", "Cooperative member") # variable names to appear in plots

coef4 <- seregfb[,1] # repeat for categories
coef4 <- unname(coef4[c(2,3)])
coef4a <- sereg5a[,1]
coef4a <- unname(coef4a[c(2,3)])
coef4b <- sereg7a[,1]
coef4b <- unname(coef4b[c(2,3)])
coef4c <- sereg9a[,1]
coef4c <- unname(coef4c[c(2,3)])
coef4e <- sereg11a[,1]
coef4e <- unname(coef4e[c(2)])
coef4g <- sereg6a[,1]
coef4g <- unname(coef4g[c(2,3)])
coef4h <- sereg16[,1]
coef4h <- unname(coef4h[c(2,3)])
coef4i <- sereg18[,1]
coef4i <- unname(coef4i[c(2,3)])
coef4j <- sereg20[,1]
coef4j <- unname(coef4j[c(2,3)])
coef4k <- sereg22[,1]
coef4k <- unname(coef4k[c(2,3)])
coef4l <- sereg24[,1]
coef4l <- unname(coef4l[c(2,3)])
coef4m <- summary(reg28)$coefficients
coef4m <- unname(coef4m[c(2,3)])
coef4n <- sereg30[,1]
coef4n <- unname(coef4n[c(2,3)])

se4 <- seregfb[,2]
se4 <- unname(se4[c(2,3)])
se4a <- sereg5a[,2]
se4a <- unname(se4a[c(2,3)])
se4b <- sereg7a[,2]
se4b <- unname(se4b[c(2,3)])
se4c <- sereg9a[,2]
se4c <- unname(se4c[c(2,3)])
se4e <- sereg11a[,2]
se4e <- unname(se4e[c(2)])
se4g <- sereg6a[,2]
se4g <- unname(se4g[c(2,3)])
se4h <- sereg16[,2]
se4h <- unname(se4h[c(2,3)])
se4i <- sereg18[,2]
se4i <- unname(se4i[c(2,3)])
se4j <- sereg20[,2]
se4j <- unname(se4j[c(2,3)])
se4k <- sereg22[,2]
se4k <- unname(se4k[c(2,3)])
se4l <- sereg24[,2]
se4l <- unname(se4l[c(2,3)])
se4m <- sqrt(diag(vcov(reg28)))
se4m <- unname(se4m[c(2,3)])
se4n <- sereg30[,2]
se4n <- unname(se4n[c(2,3)])

low4 <- coef4 - (1.645*se4)
high4 <- coef4 + (1.645*se4)
low4a <- coef4a - (1.645*se4a)
high4a <- coef4a + (1.645*se4a)
low4b <- coef4b - (1.645*se4b)
high4b <- coef4b + (1.645*se4b)
low4c <- coef4c - (1.645*se4c)
high4c <- coef4c + (1.645*se4c)
low4e <- coef4e - (1.645*se4e)
high4e <- coef4e + (1.645*se4e)
low4g <- coef4g - (1.645*se4g)
high4g <- coef4g + (1.645*se4g)
low4h <- coef4h - (1.645*se4h)
high4h <- coef4h + (1.645*se4h)
low4i <- coef4i - (1.645*se4i)
high4i <- coef4i + (1.645*se4i)
low4j <- coef4j - (1.645*se4j)
high4j <- coef4j + (1.645*se4j)
low4k <- coef4k - (1.645*se4k)
high4k <- coef4k + (1.645*se4k)
low4l <- coef4l - (1.645*se4l)
high4l <- coef4l + (1.645*se4l)
low4m <- coef4m - (1.645*se4m)
high4m <- coef4m + (1.645*se4m)
low4n <- coef4n - (1.645*se4n)
high4n <- coef4n + (1.645*se4n)

resModel4 <- data.frame(param = varnames2, 
                        low2.5 = c(low4, low4a, low4b, low4c, low4e, low4g, low4h, low4i, low4j,
                                   low4k, low4l, low4m, low4n), 
                        est = c(coef4, coef4a, coef4b, coef4c, coef4e, coef4g, coef4h, coef4i, coef4j,
                                coef4k, coef4l, coef4m, coef4n), 
                        up2.5 = c(high4, high4a, high4b, high4c, high4e, high4g, high4h, high4i, high4j,
                                  high4k, high4l, high4m, high4n)) 
resModel4$Model <- modelnames <- c("Selection", "Selection", 
                                   "No lag", "No lag", 
                                   "Drop Greece", "Drop Greece", 
                                   "Year FEs", "Year FEs", 
                                   "Disagg cooperation",  
                                   "No imputation", "No imputation",
                                   "Drop ESM", "Drop ESM",
                                   "Deference Control", "Deference Control",
                                   "Bilateral Controls", "Bilateral Controls",
                                   "Weak states", "Weak states",
                                   "World Bank control", "World Bank control",
                                   "Multilevel", "Multilevel",
                                   "Region FEs", "Region FEs")

plot4 <- ggplot(resModel4, aes(x=param, y=est, color = Model)) + 
  geom_pointrange(limits, position=position_dodge(width=0.2)) + 
  geom_hline(yintercept=0, color="red") + 
  theme_bw() + 
  theme(text = element_text(size=22), axis.title.y = element_blank(),
        axis.text.y = element_blank(), legend.position="bottom") +
  ylab(NULL) + 
  coord_flip() +
  ggtitle("Categories") +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_x_discrete(limits = rev(levels(as.factor(resModel4$param))))

g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

mylegend<-g_legend(plot4)

p3 <- grid.arrange(arrangeGrob(plot3 + theme(legend.position="none"),
                               plot4 + theme(legend.position="none"),
                               nrow=1),
                   mylegend, nrow=2, heights=c(10, 2)) # Appears as Figure 2 in paper

####################### Descriptive plots ##################################################

hist(sort(prelimimp$count), xlab = "Number of Conditions", density = 30, col = "blue",
     breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120), ylab = "Frequency", main = "",
     cex.lab=1.5, cex.axis=1.5) # histogram of conditionality to show skewness; Figure A7 in appendix

map <- as.data.frame(cbind(as.character(prelimimp$ccode), as.character(prelimimp$ctry), prelimimp$year, prelimimp$count)) # prep data for map of average conditionality burdens
colnames(map) <- c("ccode", "ctry","year", "count") # change col names
map$count <- as.numeric(as.character(map$count)) # change format
map <- aggregate(list(map$count), by = list(map$ccode, map$ctry), mean) # aggregate into mean counts
colnames(map) <- c("ccode", "ctry", "count") # change col names
maps <- joinCountryData2Map(map, joinCode = "ISO3", nameJoinColumn = "ccode", nameCountryColumn = "ctry") # assign data to map
maps <- subset(maps, continent != "Antarctica")
mapss <- mapCountryData(mapToPlot = maps, nameColumnToPlot = "count", mapRegion = "world", 
                        catMethod = "quantiles", colourPalette = brewer.pal(11, name = "Reds"),
                        addLegend = FALSE, mapTitle = "")
do.call( addMapLegend, c(mapss, legendWidth=0.5, legendMar = 5))
maps1 <- joinCountryData2Map(map, joinCode = "ISO3", nameJoinColumn = "ccode", nameCountryColumn = "ctry") # plot map
# Figure A6 in appendix

mnp <- aggregate(prelimimp$count, by = list(prelimimp$year), mean) # aggregate data into average counts over time
mnp <- mnp[mnp$Group.1 < 2015,]
plot(mnp, type = "l", xlab = "Year", ylab = "Average Number of Conditions",
     cex.lab=1.5, cex.axis=1.5) # plot average conditions over time; Figure A5 in appendix

prelimimp$change <- ifelse(prelimimp$ccode == lag(prelimimp$ccode) & 
                             prelimimp$oomem < lag(prelimimp$oomem), -1,
                           ifelse(prelimimp$ccode == lag(prelimimp$ccode) & 
                                    prelimimp$oomem > lag(prelimimp$oomem), 1, 0)) # plot to show change in oo mem
change <- prelimimp[prelimimp$change != 0,]
change$ctry <- as.character(change$ctry)
panelView(oomem~1, data = change, index = c("ctry","year"), xlab = "Year", ylab = "Country",
          ignore.treat = T, legendOff = T, main="",
          cex.lab=16, cex.axis=16) # Figure A4 in appendix

prelimimp$change <- ifelse(prelimimp$ccode == lag(prelimimp$ccode) & 
                             prelimimp$coopmem < lag(prelimimp$coopmem), -1,
                           ifelse(prelimimp$ccode == lag(prelimimp$ccode) & 
                                    prelimimp$coopmem > lag(prelimimp$coopmem), 1, 0)) # plot to show change in coop mem
change <- prelimimp[prelimimp$change != 0,]
change <- change[-c(1,4,7),] # delete NA and duplicate rows
change$ctry <- as.character(change$ctry)
panelView(coopmem~1, data = change, index = c("ctry","year"), xlab = "Year", ylab = "Country",
          ignore.treat = T, legendOff = T, main="",
          cex.lab=16, cex.axis=16) # Figure A3 in appendix

rfash <- read.csv("rfalendshare.csv") # load data on RFA lending, IMF lending to RFA members, and 3 year moving averages
rfashp <- read.csv("rfasharep.csv") # load data with just the moving average shares manually copied from prior sheet

####################### Descriptive stats ########################
covs_desc <- c("Duration", "Number of conditions", "PRGF", "GDPPC", "U.S. aid", "Openness", 
               "Debt service / exports", "Short-term debt / exports", "FDI / GDP", "Quota",
               "UN voting (ideal pt dist from U.S.)", "Polity2", "Checks", "Time from last IMF program",
               "Liberal ideology", "Cooperative member", "Outside option member", 
               "Cooperative amount", "Outside option amount", "Election year",
               "War", "UNSC member", "Number of categories") # covariate labels
desc <- read.csv("prelimimp.csv")
desc$stone <- rowSums(desc[,c("DEB","FIN","FP","EXT","RTP","SOE","LAB","PRI","SP","POV",
                                        "INS","ENV","OTH")]>0) # create number of categories variable
desc <- desc[,-c(1:6, 8:9, 11:24, 26, 29, 33:34, 36, 39, 44:53, 58:59)]
stargazer(desc, omit.summary.stat = c("p25", "p75"), digits = 2, style = "ajps",
          covariate.labels = covs_desc) # Table A10 in appendix

prelim <- prelim[!duplicated(prelim$ID),]
prelim$loansize <- gsub("(.*),.*", "\\1", prelim$loansize) # remove values after commma
prelim$loansize <- as.numeric(as.character(prelim$loansize)) # change format
prelim$loansize[is.na(prelim$loansize)] <- 0 # make NA loan size zero to allow aggregation
prelimagg <- aggregate(prelim$loansize, by = list(prelim$ctry, prelim$year), FUN = sum) # sum IMF loans by country-year
colnames(prelimagg) <- c("ctry", "year", "imfloans") # change column names
prelimaggflar <- prelimagg[prelimagg$ctry %in% c("Bolivia", "Colombia",	"Costa Rica",	"Ecuador",	
                                                 "Paraguay",	"Peru",	"Uruguay",	"Venezuela"),] # restrict to FLAR countries
prelimaggflar <- aggregate(prelimaggflar$imfloans, by = list(prelimaggflar$year), FUN = sum) # aggregate by year 
prelimaggamf <- prelimagg[prelimagg$ctry %in% c("Jordan",	"Emirates",	"Bahrain",	"Tunisia",	"Algeria",	
                                                 "Djibouti",	"Saudi Arabia",	"Sudan",	"Syria",	
                                                 "Somalia",	"Iraq",	"Oman",	"Qatar",	"Comoros",	
                                                 "Kuwait",	"Lebanon",	"Libya",	"Egypt",	"Morocco",	
                                                 "Mauritania",	"Yemen"),] # restrict to AMF countries
prelimaggamf <- aggregate(prelimaggamf$imfloans, by = list(prelimaggamf$year), FUN = sum) # aggregate by year
prelimaggefsd <- prelimagg[prelimagg$ctry %in% c("Russia", "Armenia", "Belarus", "Kazakhstan",
                                                 "Kyrgyzstan", "Tajikistan"),] # restrict to EFSD countries
prelimaggefsd <- aggregate(prelimaggefsd$imfloans, by = list(prelimaggefsd$year), FUN = sum) # aggregate by year
prelimaggbop <- prelimagg[prelimagg$ctry %in% c("Bulgaria",	"Croatia",	"Czech Republic",	"Denmark",	
                                                "Hungary",	"Latvia",	"Poland",	"Romania",	"Sweden",	
                                                "United Kingdom"),] # restrict to BoP countries
prelimaggbop <- aggregate(prelimaggbop$imfloans, by = list(prelimaggbop$year), FUN = sum) # aggregate by year
prelimaggesm <- prelimagg[prelimagg$ctry %in% c("Austria", "Belgium", "Cyprus", "Estonia", "Finland", 
                                                "France", "Germany", "Greece", "Ireland", "Italy", 
                                                "Latvia", "Lithuania", "Luxembourg", "Malta", 
                                                "Netherlands", "Portugal", "Slovak Republic", 
                                                "Slovenia", "Spain"),] # restrict to ESM countries
prelimaggesm <- aggregate(prelimaggesm$imfloans, by = list(prelimaggesm$year), FUN = sum) # aggregate by year
prelimaggefsm <- prelimagg[prelimagg$ctry %in% c("Austria", "Belgium", "Bulgaria", "Croatia", "Cyprus", 
                                                 "Czech Republic", "Denmark", "Estonia", "Finland", 
                                                "France", "Germany", "Greece", "Hungary", "Ireland", "Italy", 
                                                "Latvia", "Lithuania", "Luxembourg", "Malta", 
                                                "Netherlands", "Poland", "Portugal", "Romania", "Slovak Republic", 
                                                "Slovenia", "Spain", "Sweden", "United Kingdom"),] # restrict to EFSM countries
prelimaggefsm <- aggregate(prelimaggefsm$imfloans, by = list(prelimaggefsm$year), FUN = sum) # aggregate by year
# The data above appear in Tables A4--A9 and were manually inputted

####################### Text analysis ########
##->## Cleaning and preprocessing
conditions <- read.csv("IMF_cond.csv") # load conditions tab from Kentikelenis et al. (2016)
Corpus <- Corpus(VectorSource(conditions$Text)) # convert to corpus
Corpus.clean <- tm_map(Corpus, stripWhitespace) # strip white space
Corpus.clean <- tm_map(Corpus.clean, removePunctuation) # remove punctuation
Corpus.clean <- tm_map(Corpus.clean, tolower) # convert to lower case
dataframe <- data.frame(text=unlist(sapply(Corpus.clean, `[`)), stringsAsFactors=F) # convert to frame
conditions <- cbind(conditions, dataframe) # re-append text to frame
imfsub <- read.csv("prelimimp.csv") # load data on oo and coop IO membership
imfsub <- imfsub[c(6, 54, 55)] # subset to just ID and membership columns
colnames(conditions)[4] <- "ID" # harmonize ID column name with other frame
conditions <- merge(unique(conditions), unique(imfsub), by = c("ID"), all.x = T) # merge cooperation data into text data
coop <- conditions[conditions$coopmem==1,] # frame with just coop observations
comp <- conditions[conditions$oomem==1,] # frame with just comp observations
processed <- textProcessor(conditions$text, metadata = conditions) # process text data
out <- prepDocuments(processed$documents, processed$vocab, processed$meta) # build data

##->## Topic modeling for oo member treatment
imf.multi <- manyTopics(processed$documents, processed$vocab, K = 6:16, 
                        prevalence = ~ coopmem + oomem, max.em.its = 100, runs = 10, 
                        data = processed$meta, seed = 1234) # determine how many topics to use
t <- imf.multi$semcoh # prepare data to plot
s <- imf.multi$exclusivity
sc <- c()
ex <- c()
cols <- c()
for(i in 1:11){
  sc <- c(sc, mean(t[[i]]))
  ex <- c(ex, mean(s[[i]]))
}
plot(sc, ex, xlim = c(-100, -205),
     ylim = c(8.0, 10.0), color = cols) # the plot suggests that around 10 topics is optimal

imf.model <- selectModel(processed$documents, processed$vocab, K=10, 
                          prevalence = ~ oomem + coopmem, max.em.its = 100,
                          runs=10, data = processed$meta, seed = 1234) # run model with 10 topics for becoming oomem or coopmem
plotModels(imf.model) # 2 is better
imf.model1 <- imf.model$runout[[2]] # select 2
labelTopics(imf.model1) # see most common words for labeling
imf.model1.effect <-  estimateEffect(1:10 ~ oomem + coopmem, imf.model1, 
                                     meta=processed$meta, uncertainty = "Global")
summary(imf.model1.effect) # examine which topics are significant

plot(imf.model1.effect, covariate = "oomem", topics = c(1,2,5,7), 
     method = "difference", cov.value1=1, cov.value2=0, ci.level = 0.95, 
     xlab = "Change in topical prevalence when joining competitive IO", 
     main = "", xlim = c(-.05,.05), labeltype = "custom", 
     custom.labels = c('Tax law','Publications',
                       'Privatization','Bureaucratic tasks'),
     cex.lab=1.5, cex.axis=1.5) # plot topical prevalence for oo member treatment; Figure 3 in text

oo_1 <- findThoughts(imf.model1, texts=processed$meta$text, topics=1, n=15) # look at which topics are representative for each significant topic; the # sign demarcates the ones that will appear in the paper
oo_1$index # Tax law (-)
processed$meta$text[10269] # no granting of exemptions to any company [...] and no renewal of existing exemptions
processed$meta$text[11536] # adopt the new income tax code
processed$meta$text[2818] # adopt [...] tax reform package including elimination of tax exemptions and preferential regimes

oo_2 <- findThoughts(imf.model1, texts=processed$meta$text, topics=2, n=15) # look at which topics are representative for each significant topic; the # sign demarcates the ones that will appear in the paper
oo_2$index # Publications (+)
processed$meta$text[8480] # sign memoranda of understandings to delineate responsibilities
processed$meta$text[1200] # publish financial statements on website

oo_5 <- findThoughts(imf.model1, texts=processed$meta$text, topics=5, n=15) # look at which topics are representative for each significant topic; the # sign demarcates the ones that will appear in the paper
oo_5$index # Privatization (-)
processed$meta$text[7034] # conclusion of at least six contracts for the privatization of either large companies or pools of six companies
processed$meta$text[30350] # privatization of the drilling company sonafor the tapestry weaving company msad and the reinsurance company senre

##->## Topic modeling for coop member treatment
plot(imf.model1.effect, covariate = "coopmem", topics = c(1,2,3,4,7,8,9), 
     method = "difference", cov.value1=1, cov.value2=0, ci.level = 0.95, 
     xlab = "Change in topical prevalence when joining cooperative IO", 
     main = "", xlim = c(-.15,.15), labeltype = "custom", 
     custom.labels = c('Tax law', 'Publications', 'Public spending', 'External debt payments',
                       'Foreign exchange policy', 'Bureaucratic tasks',
                       'Trade policy')) # plot topical prevalence for coop member treatment; Figure 4 in text

plot(imf.model1, type="labels", n = 10, topics=c(1:5,7:10), topic.names = c('Tax law',
                                                                            'Publications',
                                                                            'Public spending',
                                                                            'External debt payments',
                                                                            'Privatization',
                                                                            'Foreign exchange policy',
                                                                            'Bureaucratic tasks',
                                                                            'Trade policy',
                                                                            'Monetary policy'
                                                                            )) # plot key words from conditions (Table A10)

####################### IMF Network Diagram #####

nodes <- read.csv("IMF-NODES.csv", header=T, as.is=T) # hand-coded data on collaboration amenable to network analysis
links <- read.csv("IMF-EDGES.csv", header=T, as.is=T)
net <- graph_from_data_frame(d=links, vertices=nodes, directed=F) 

# Generate colors based on cooperation type:
colrs <- c("gray50", "tomato", "gold")
V(net)$color <- colrs[V(net)$type]

# Compute node degrees (#links) and use that to set node size:
deg <- igraph::degree(net, mode="all")
V(net)$size <- deg

# The labels are currently node IDs.
# Setting them to NA will render no labels:
V(net)$label <- NA

# Set edge width based on weight:
E(net)$width <- E(net)$weight

#change arrow size and edge color:
E(net)$arrow.size <- .2
E(net)$edge.color <- "gray80"

# We can even set the network layout:
graph_attr(net, "layout") <- layout_with_lgl
plot(net, vertex.label=V(net)$media, layout = layout_in_circle) 

# Plot different types of cooperation separately
net.m <- net - E(net)[E(net)$type=="cofinance" | E(net)$type=="other"] # another way to delete edges:
net.h <- net - E(net)[E(net)$type=="information" | E(net)$type=="other"] # using the minus operator

# Plot the two links separately:
par(mfrow=c(1,2))
plot(net.h, vertex.color="orange", layout=layout_in_circle, main="Co-financing", vertex.label=V(net)$media)
plot(net.m, vertex.color="lightsteelblue2", layout=layout_in_circle, main="Information Sharing", vertex.label=V(net)$media)
# Figure 1 in paper

## Plot cooperation over time
par(mfrow=c(1,1))
imfcollab <- read.csv("imf_collab.csv")
mnc <- aggregate(imfcollab$coopagg, by = list(imfcollab$year), sum) # aggregate data into average counts over time
plot(mnc, type = "l", xlab = "Year", ylab = "Total Cooperation",
  cex.lab=1.5, cex.axis=1.5) # Figure A1 in appendix

## Plot competition over time
par(mfrow=c(1,1))
mno <- aggregate(prelimimp$oomem, by = list(prelimimp$year), sum) # aggregate data into average counts over time
plot(mno, type = "l", xlab = "Year", ylab = "Total Competition",
     cex.lab=1.5, cex.axis=1.5) # Figure A2 in appendix


####################### Which observations are missing? ####
miss <- prelimlag[c(2,3,4,10,25,27,28,30,31,32,35,38,40,41,42,43,54:57,60:63)] # create dataset just with covariates used for analysis
miss$complete <- complete.cases(miss) # create variable for whether observation is complete i.e. no missingness
colMeans(is.na(miss)) # calculate share missing by variable. I should leave proportion neoliberal out of the no imputation model at minimum
miss <- miss[-c(8,9,15)] # drop debt service / exports, short-term debt / exports, and prop neoliberal
miss$complete <- complete.cases(miss) # re-create variable for whether observation is complete i.e. no missingness
miss$incomplete <- 1-miss$complete # recode as incomplete
missctry <- aggregate(miss$incomplete, list(miss$ctry, miss$ccode), FUN = sum) # calculate number of obs missing by country
colnames(missctry) <- c("ctry", "ccode", "numincomplete") # rename cols
missyr <- aggregate(miss$incomplete, list(miss$year), FUN = sum) # calculate number of obs missing by year
colnames(missyr) <- c("year", "numincomplete") # rename cols
missyr <- missyr[-c(38:39),] # cut 2015-16, which are not in the data set

maps <- joinCountryData2Map(missctry, joinCode = "ISO3", nameJoinColumn = "ccode", nameCountryColumn = "ctry") # assign data to map
maps <- subset(maps, continent != "Antarctica")
mapss <- mapCountryData(mapToPlot = maps, nameColumnToPlot = "numincomplete", mapRegion = "world", 
                        catMethod = "quantiles", colourPalette = brewer.pal(5, name = "Reds"),
                        addLegend = FALSE, mapTitle = "")
do.call( addMapLegend, c(mapss, legendWidth=0.5, legendMar = 5)) # Figure A8 in appendix

plot(missyr, type = "l", xlab = "Year", ylab = "Number of Missing Observations",
     cex.lab=1.5, cex.axis=1.5) # plot missing observations over time; Figure A9 in appendix

####################### Network regression ####
## Create datasets
imfcollab <- read.csv("imf_collab.csv") # load hand-coded data on collaboration among emergency lending IOs
imfcollab$coopagg <- imfcollab$cofin + imfcollab$infoshare # code aggregate cooperation to only include inforshare and cofin
imfcollab$coopbin <- 0 # create a binary for cooperation
imfcollab$coopbin[imfcollab$coopagg > 0] <- 1 # code the binary
imfcollab2 <- read.csv("imf_collab_clean.csv") # load data with hand-coded covariates
imfcollab2$memdiff <- abs(imfcollab2$io1nummem - imfcollab2$io2nummem) # calculate differences in memberships
imfcollab2$dyid <- paste(imfcollab2$io1short, imfcollab2$io2short) # assign dyad ID
lagdv1 <- slide(imfcollab2, Var = "infoshare", GroupVar = "dyid",
                slideBy = 1, NewVar = 'infosharelag') # create lag variable for infoshare
lagdv2 <- slide(imfcollab2, Var = "cofin", GroupVar = "dyid",
                slideBy = 1, NewVar = 'cofinlag') # create lag variable for cofin
imfcollab2 <- cbind(imfcollab2, lagdv1$infosharelag) # merge in first lag
imfcollab2 <- cbind(imfcollab2, lagdv2$cofinlag) # merge in second lag
colnames(imfcollab2)[25] <- "infosharelag" # rename columns
colnames(imfcollab2)[26] <- "cofinlag"
un <- read.csv("unidealpt.csv") # load dyadic UN voting data downloaded from Bailey et al. (2017) in November 2018 V21 https://doi.org/10.7910/DVN/LEJUQZ
imfcollab2$imfbin <- ifelse(imfcollab2$io1short == "IMF", 1, 0) # create binary variable for World Bank involvement
imfcollab2$io1topsh <- as.character(imfcollab2$io1topsh) # change data format to character
imfcollab2$io2topsh <- as.character(imfcollab2$io2topsh)
imfcollab2$hqctry1 <- as.character(imfcollab2$hqctry1)
imfcollab2$hqctry2 <- as.character(imfcollab2$hqctry2)
imfcollab2$io1topshc <- countrycode(imfcollab2$io1topsh, origin = "country.name", destination = "cown") # create country code variables for both shareholder and HQ countries 
imfcollab2$io2topshc <- countrycode(imfcollab2$io2topsh, origin = "country.name", destination = "cown")
imfcollab2$idealptsh <- NA # create empty variables for UN distance, alliance ties, and rivalry-peace between leading stakeholders
imfcollab2$idealpthq <- NA # variables are then generated in the loops below
for (i in 1:181) {
  un1 <- un$idealpt[un$ccode == imfcollab2$io1topshc[i] & un$year == imfcollab2$year[i]]
  un2 <- un$idealpt[un$ccode == imfcollab2$io2topshc[i] & un$year == imfcollab2$year[i]]
  imfcollab2$idealptsh[i] <- ifelse(is.numeric(abs(un1 - un2)), abs(un1 - un2),NA)
  imfcollab2$idealptsh[i] <- ifelse(imfcollab2$io1topshc == imfcollab2$io2topshc, 0, imfcollab2$idealptsh[i])
} # loop to complete UN voting variables

## Analysis
r1 <- glm.nb(infoshare ~ idealptsh + factor(year), data = imfcollab2)
(ser1 <- coeftest(r1, vcov = vcovHC(r1, type = "HC1", cluster = imfcollab2$dyid))) # calculate SEs
r2 <- glm.nb(cofin ~ idealptsh + factor(year), data = imfcollab2)
(ser2 <- coeftest(r2, vcov = vcovHC(r2, type = "HC1", cluster = imfcollab2$dyid))) # calculate SEs
stargazer(r1,r2, se = list(ser1[,2], ser2[,2]),
          dep.var.labels = c("Information Sharing", "Co-financing"),
          covariate.labels = c("UN voting (ideal pt dist)", "Year"),
          style = "ajps", omit = c("Constant","year"), add.lines = yearfe,
          omit.stat = c("ll", "aic", "theta")) # Table A11 in appendix

sink()

